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We present an extensive first-principles study of the pressure dependence of the formation enthalpies of all 
the know vacancy and self-interstitial configurations in silicon, in each charge state from —2 through +2. The 
neutral vacancy is found to have a formation volume that varies markedly with pressure, leading to a remarkably 
large negative value (—0.68 atomic volumes) for the zero-pressure formation volume of a Frenkel pair (V + 1). 
The interaction of volume and charge was examined, leading to pressure-Fermi level stability diagrams of the 
defects. Finally, we quantify the anisotropic nature of the lattice relaxation around the neutral defects. 



PACS numbers: 61.72.Bb,61.72Ji,66.30.Hs 

INTRODUCTION 

Nearly perfect crystals of silicon are of great technological 
importance, yet silicon self-diffusion is still not completely 
understood. Unlike the situation in metals, the equilibrium 
concentrations of vacancies and self-interstitials in Si are be- 
lieved to comparable, and very low, making detection of them 
problematic. Experimental data is fragmentary, and simula- 
tions do not all agree. Controversy remains over the relative 
importance of vacancies and interstitials to self-diffusion at 
different temperatures and the relative magnitudes of the mi- 
gration enthalpy to the formation enthalpy of these defects 
[1,2]. We set aside the possibility of diffusion by a direct ex- 
change mechanism [3] due to the low prefactor that has been 
calculated [4]. Then the self-diffusivity of silicon is the sum 
of the diffusion of Si due to vacancies and due to interstitials, 

D Si =c v D v + c I D I (1) 

where the atomic fraction cx = Cx/Csi and Csi = 
5.00 x 10 22 cm~ 3 . The contribution of vacancies to Si dif- 
fusion is proportional to the concentration of vacancies Cy 
and the diffusivity of vacancies Dy, and likewise for inter- 
stitials. The equilibrium concentration of a defect X is C x q = 
Csi exp (—gx/ksT), where g x is the Gibbs free energy of for- 
mation of one defect. The diffusivity of a vacancy can be 
written as Dy = (£/6)A 2 Voexp(— gy/ksT), where £ is the 
coordination number, A is the bond length, Vo is an attempt 
frequency, and gy is the Gibbs free energy for the vacancy to 
exchange with one of its neighbors. The diffusivity of self- 
interstitials can be written similarly, but with a different geo- 
metric factor. Of course, G = H — TS, and if the entropy and 
enthalpy are assumed to be constant with respect to tempera- 
ture, entropy may be combined with the pre-exponential fac- 
tor, leaving only the enthalpy as a model parameter in the ex- 
ponent. 



Recent isotope tracer experiments [2] fit equilibrium silicon 
self-diffusivity to a single Arrhenius term 

with D c s f = 530+170 cm 2 /s and h Si = 4.75 ± 0.04 eV, sug- 
gesting that either vacancies or self-interstitials dominate self- 
diffusion over the entire temperature range studied — or in- 
stead that they switch over from one to the other but with sim- 
ilar values of h^ . Of course, vacancies and self-interstitials 
diffuse by exchanging with lattice atoms, and thus cannot be 
isotopically tagged. 

The concentrations C and diffusivities D of vacancies in sil- 
icon are difficult to measure separately with any accuracy, and 
likewise with self-interstitials. Estimated equilibrium trans- 
port capacities c^Dx = d x exp (— (li x + h x ) /IcbT) have been 
derived from experimental studies of metal diffusion in sili- 
con, the most recent of which report values of h v + hy ranging 
from 4.03 eV to 4. 14 eV, and h\ + hf from 4.84 eV to 4.95 eV 
[2, 5]. On the other hand, the latest published work utiliz- 
ing dopant diffusion arrives at an estimate of h d = h v + hy = 
4.86 eV and h d = h\ + hf = 4.68 eV [1], contrary to the long- 
held assumption that hf > hy. Ion implantation, thermal ox- 
idation, and nitridation increase the concentration of intrinsic 
defects above their equilibrium concentrations (Cx 3> C x q ), in- 
creasing the self- and dopant diffusivity; however, measuring 
activation enthalpies from such experiments requires assump- 
tions about traps and other simplifications to differentiate be- 
tween interstitial and vacancy mechanisms. 

First-principles methods can be used to separately calculate 
formation and migration enthalpies, among other quantities. 
The enthalpies can be obtained from the energies calculated 
at different volumes as H = E + PV + qe?, where £f is the 
Fermi level and q is the charge of the defect in electron units. 
Since V = dH/dP, the formation volume v x tells us how the 



formation enthalpy h x , and thus the equilibrium concentra- 
tion varies with pressure. Similarly, v™ tells us how much 
pressure enhances or retards the migration of a defect once it 
has entered the lattice. This is of technological interest be- 
cause of the large stresses and strains that exist near the sur- 
faces and interfaces of silicon-based integrated circuits (e.g. 
heteroepitaxial growth of Si upon Sii-jc-vGe^Cy or dielectric 
substrates), which may have a strong effect on concentrations 
and diffusion of intrinsic defects. 

METHODS 

Total energies were calculated using density functional the- 
ory (DFT), as implemented in the code VASP [6-9]. All calcu- 
lations were performed with the PW91 exchange-correlation 
functional. Ion cores were represented with Vanderbilt ultra- 
soft pseudopotentials, allowing plane-wave energy cutoffs of 
HRy = 150 eV. The Brillouin zone was sampled with k 
points equivalent to a 4x4x4 Monkhorst-Pack mesh in a con- 
ventional cubic (8-atom) cell. Periodic boundary conditions 
were used with primitive (2-atom) silicon cells for calculation 
of pure silicon and larger (mostly 128- or 256-site) supercells 
for calculations involving defects. 

Our reference point was perfect silicon in the diamond cu- 
bic structure. The total energies at different volumes were fit 
to a Birch-Murnaghan equation of state [10]. In our calcu- 
lations, silicon had a cohesive energy of 4.53 eV and equi- 

° 3 

librium atomic volume vs; = 20.34 A . The bulk modulus 
at P = was Bq = 0.88 Mbar and its pressure derivative was 
B;=4.02. 

Similar calculations were performed with supercells con- 
taining defects. In each case, the supercells were set at a par- 
ticular volume (scaled isotropically from the perfect lattice) 
and the ionic coordinates were fully relaxed to build up a list 
of at least seven energy-volume data points. Then the same 
equation of state was fit to the energy-volume data to find 
the enthalpy-pressure relationship. We define as the to- 
tal enthalpy and vs; as the volume of a silicon atom in the 
perfect crystal. A supercell containing 256 lattice sites and 
a vacancy includes 255 silicon atoms, so we define the for- 
mation enthalpy of a vacancy as h v = i/y[Si25s] — 255/isi, 
and likewise the formation volume is v v = Vy [Si25s] — 255 vs; . 
Similarly, for self-interstitials = Hj[Si25i] — 251 and 
vjr = Vy[Si257] — 257i'si . On the other hand, the quanti- 
ties of interest in elasticity are the relaxation volumes v™ 1 = 
Vx[Si256±i] — 256 vs;. The defects described in this work are 
shown in Figure 1 and will be commented on later. 

VACANCIES 

Intuitively, the simplest point defect is a neutral lattice va- 
cancy, which we will label V®: removing an atom from a per- 
fect lattice. However, the neighboring atoms will tend to re- 
bond in ways that make the defect less symmetric, particu- 



2 




(d) V B (e) I T (f) I H 



FIG. 1 : Geometries of defects in relation to a conventional cubic 
unit cell of Si 



a (A) 

5.55 5.5 545 54 5.35 




-60 -40 -20 20 40 60 



P (kbar) 

FIG. 2: Jahn-Teller splitting of energy eigenvalues. Shown are the 
band edges and the defect levels for V® with (gray) and V® with 
Dod (black) symmetry. Pressures are approximate. 

larly in a covalently bonded crystal like silicon. If the atoms 
are forced to maintain a Tj symmetry, the four neighbors draw 
in toward the center, pulling the rest of the lattice with them. 
However, the ground state involves a Tj — > Z^/ symmetry- 
breaking relaxation explained in the early days of quantum 
chemistry by Jahn and Teller [11]. The Jahn-Teller distortion 
of the neutral lattice vacancy in Si is now well-established by 
experiment [12] and theory [13]. We can see in Figure 2 the 
splitting of the triply degenerate Ti level (occupied by two 
electrons) into a filled lower level and two degenerate empty 
upper levels, all still in the band gap. 

Figure 3 shows that the Jahn-Teller distortion reduces the 
enthalpy of the defect by about 0.25 eV at P = 0, and results 
in a further contraction of 0.40vs;. As a result, the formation 
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TABLE I: Vacancy formation enthalpies and volumes at P = and 



£p — £{. 
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h f (eV ) 


AvsO 


h s (eV ) 


v f (v Si ) 


-2 


4.33 


0.08 


4.14 


-0.27 


-1 


3.87 


0.01 


3.89 


-0.30 





3.69 


-0.07 


3.97 


-0.30 


+ 1 
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0.21 


4.29 


-0.37 


+2 


4.55 


0.43 


4.90 


-0.42 



FIG. 3: Effect of pressure on formation enthalpy of neutral va- 
cancies. Shown here are curves for V® with both Tj symmetry (gray) 
and Dtj {heavy), as well as V^j (dashed). Fermi level £f fixed at in- 
trinsic level £j. 
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FIG. 4: Effect of supercell size on neutral vacancy (V°) formation 
volume at P = and £p = £j. 



volume of a T ( ) vacancy is positive in our calculations, but that 
of a D2d vacancy is a small negative value, — 0.07vsi at P = 0, 
and vanishes at about —6 kbar. This surprising result implies 
that introducing vacancies into silicon (at P = 0) actually re- 
duces the volume of the system and increases its density. A 
lattice with no relaxation would have vt = +Vsi, while per- 
fect relaxation (as with an incompressible liquid) would give 

V f y = 0. 

The Jahn-Teller splitting vanishes under sufficiently large 
tensile strain, but is approximately constant under compres- 
sive strain. Other first-principles calculations of V® (though 
performed only for P = 0) have also reported small negative 
formation volumes [14]. In simulations reported in the litera- 
ture, this distortion is stable only when using supercells with 
more than 128 atoms [15]. Our own calculations confirm this, 
with a transition from tetrahedral to tetragonal geometry oc- 
curring between 128 and 216 atoms (Figure 4). 

Another high-symmetry configuration of N — 1 atoms in a 
crystal with N sites is a so-called split vacancy, where one 
atom is at the bond center between two empty sites, a config- 
uration we label Vb- This can easily be seen as a transition 
point for vacancy migration. Table I shows that Vb has a fairly 
constant large negative formation volume, about — 0.30vs;. 



However, in contrast to the lattice vacancy, this is not due to 
a Jahn-Teller symmetry breaking. As a result, the formation 
volume of the six-coordinated Vg shows little change with su- 
percell size (Figure 4) or pressure (Figure 3). From these ar- 
guments it also follows that the vacancy migration enthalpy is 
fairly constant for positive pressure, but increases for negative 
pressures. 

We have also calculated the formation enthalpies and vol- 
umes of charged vacancies (Table I), since charged defects 
are expected to play an important role in doped silicon. Fig- 
ure 5(a) shows that the neutral lattice-centered vacancy V® 
is the stable configuration over the broadest range of Fermi 
levels within the band gap (including the intrinsic level £; = 
(e v + £ c )/2) and has a formation enthalpy of 3.69 eV at P = 0. 
The bond-centered split vacancy Vg is 0.27 eV higher. 

As mentioned above, Vb is the transition point for a vacancy 
to migrate from one Vl configuration to the next. However, 
our calculations reveal a surprising twist: The negative split 
vacancy Vg has almost the same energy as the negative lat- 
tice vacancy V£ , and at high Fermi levels, the ground state is 
in fact the split vacancy (the doubly negative split vacancy is 
even lower in energy than the lattice vacancy). This phenome- 
non was predicted decades ago by Bourgoin and Corbett [16], 
but is not well-known by all who work with silicon diffusion, 
since most industrially useful processes involve elevated tem- 
peratures. This reversal may be due to the greater number of 
bonds that can accept extra electrons in the case of the split va- 
cancy. Since the displaced atom in the split vacancy has four 
valence electrons and six neighbors, giving the system two 
extra electrons allows this atom to form six equally strained 
bonds with all of its neighbors. In other words, the minimum- 
enthalpy geometry in one charge state (—2) is a saddle point 
in a different charge state (0), and vice versa. This crossover 
in potential energy surface is the requirement for the Bourgoin 
mechanism of athermal (electronically or optically activated) 
diffusion to take place [16, 17]. Experiments have demon- 
strated that vacancies in n-type silicon can diffuse at room 
temperature or even cryogenic temperature when subjected to 
optical or electronic excitation [12, 18]. 
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FIG. 5: Intrinsic defect formation enthalpies vs. Fermi level at P — 0. The heavy solid line indicates the lowest formation enthalpy h f for a 
defect at that Fermi level. The heavy dashed line indicates the second-lowest enthalpy, which is an estimate of h d , the activation enthalpy of 
the transport capacity. The difference between the lines would then be the activation enthalpy of migration, h m . Note that slope is proportional 
to the charge of the defect. A vertical line marks the calculated location of the intrinsic Fermi level £;. 



SELF-INTERSTITIALS 

We have calculated the enthalpies of the three interstitial ge- 
ometries found to have the lowest enthalpies in previous work: 
tetrahedral (It), hexagonal (/#), and split- (110) (/(no))- Our 
self-interstitial calculations are summarized in Table II. In 
contrast to the strong pressure and Fermi level sensitivity of 
vacancies, self-interstitials present a somewhat simpler pic- 
ture. The lowest enthalpy interstitial at the intrinsic Fermi 
level and P = is Ij + , with a formation enthalpy of 3.68 eV 
and a formation volume of — 0.51vs;. (See the Methods sec- 
tion for an explanation of how formation enthalpies and vol- 
umes are defined in our calculations.) 

We can explain the stability region of the various self- 
interstitial configurations, by plotting their formation en- 
thalpies as a function of Fermi level (Figure 5(b)). We see that 
the Ij + configuration is stable mainly in /5-type Si, while at 
slightly elevated Fermi levels (« doping) the fL split inter- 
stitial is stabilized. A strong negative-C/ effect is seen, where 
the interstitial is nowhere stable in the +1 charge state. 

Much like the situation with vacancies, high-symmetry in- 
terstitial geometries can be saddle points for interstitial mi- 
gration. In particular, the hexagonal interstitial is usually as- 
sumed to be the migration point for interstitial diffusion. Fig- 
ure 5(b) suggests that this is true for Si at almost any dop- 

o 

(110) —"-H 

p- 

leads 



ing; hence interstitial migration in «-doped Si (if, 1fl \ — ► /? 
/?,,„,) has a migration barrier of only 0.06 eV, while in 



doped Si the migration mechanism Ij. 



to a barrier that increases linearly with the doping level, thus 
causing a drastic slow-down of interstitial migration. 

A particularly interesting regime is that of the Fermi level 



TABLE II: Self-interstitial formation enthalpies and volumes at 

P = and e F = £;. (vf = -v si + vf 1 ) 







v f (v Si ) 




v f (v Si ) 


I 9 


Vsi) 


1 


h f (eV ) 


h f (eV ) 


h f (eV ) \ 


-2 


5.37 


-0.73 


5.12 


-0.49 


4.70 


-0.45 


-1 


4.61 


-0.67 


4.34 


-0.43 


4.06 


-0.42 





4.06 


-0.63 


3.79 


-0.38 


3.73 


-0.41 


+ 1 


3.73 


-0.60 


4.15 


-0.41 


4.19 


-0.45 


+2 


3.68 


-0.51 


4.85 


-0.49 


4.92 


-0.51 



being close to the midgap, where the Ij + is only slightly lower 
in energy than the neutral 7° . It is easy to see that under 
such conditions, addition of electrons e.g.through electron ir- 
radiation, can lead to a mechanism without any activation bar- 
rier. This so-called Bourgoin mechanism of electrically, rather 
than thermally, activated diffusion is observed experimentally 
at cryogenic temperatures during electron irradiation of intrin- 
sic Si. One of the implications of this is that annealing after 
ion implantation should be significantly faster in n-type sili- 
con than in p-type. 



FRENKEL PAIRS 

Combining the formation volumes of the most stable 
vacancy and self-interstitial configurations indicates that a 
Frenkel pair (V + 1) should have a formation enthalpy of 
7.39 eV and a formation volume of — 0.68vsi. At elevated 
Fermi levels (i.e. high n doping), the stable self-interstitial 
species becomes P? no \> which would lead to a Frenkel pair 
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formation enthalpy of 7.44 eV and formation volume of 
-0.48v si . 

On the other hand, Huang diffuse x-ray scattering experi- 
ments performed by Ehrhart et al. [19, 20] have been taken 
to conclude that the formation volumes of Frenkel pairs (pro- 
duced by 2.5 MeV e~ or 4.5 keV He implantation) are quite 
small, 0.1 vsi or l ess - This result is based on the assumption 
that at cryogenic temperatures, electron irradiation generates a 
substantial number of Frenkel pairs with a separation distance 
of only about 8 A. 

The disagreement of the Erhart experiments with our first- 
principles calculations suggests that Frenkel pairs with short 
separation distances behave differently from the isolated de- 
fects. This is not an unreasonable theory when considering 
our previous findings (see section Vacancies) of the sensitiv- 
ity of the formation volume of the vacancy to its concentra- 
tion (see Fig. 4). To see whether vacancies and interstitials 
at close proximity can behave anomalously, we constructed 
supercells containing a self-interstitial and a vacancy approx- 
imately 8 A apart. Two different interstitial geometries were 
used, It and /(no)- Only neutral supercells were considered, 
since the claim of Ehrhart et al. is that irradiation of silicon 
creates large concentrations of Frenkel pairs that are either en- 
tirely neutral or form donor-acceptor pairs. In this case, the 
Frenkel pair would not consist of V® but perhaps of 

+-f? 110 \> V£ +^7' or VL +^t :+ - Formation volumes are 
calculated as before. The results are shown in Table III, to- 
gether with summed values for the isolated defects. 



TABLE III: Frenkel pair enthalpies and volumes of formation. 

Values for isolated defects are summed for convenience. 



geometry 


r(A) 


//(eV) 


v f (v Si ) 


(VbIt) 


5.88 


6.88 


-0.89 


(VlIt) 


8.83 


7.18 


-0.63 


(v L i im )° 


8.21 


7.31 


-0.48 


V°+/++ 
v L -+/++ 


oo 
oo 
oo 


7.44 
7.39 
8.04 


-0.48 
-0.68 
-0.62 



The formation volume for V£ +/^ 110 ) i s nearly the same 
whether both defects are in the same supercell or calculated 
separately. The formation enthalpy is a bit lower for the pair 
than for the isolated defects, which indicates the strength of 
the their attraction. Interestingly, in the case of iy L I T ) { \ the 
relaxation of the atoms surrounding the vacancy took the form 
of a bent square, rather than two pairs seen in the isolated va- 
cancy, even when the symmetry was broken in the direction of 
the expected geometry. In summary, for either the T or (1 10) 
geometry, the Frenkel pair formation volume is a sizable neg- 
ative number, contrary to the experimental results of Erhart et 
al. 

In light of this, the assertion of Erhart et al. that irradia- 
tion of Si at cryogenic temperatures will create a large num- 



ber of Frenkel pairs with short separation distances is incon- 
sistent with first-principles calculations. Instead, our calcula- 
tions suggest that far more defect clustering occurs than would 
be expected from purely classical migration mechanism of 
neutral defects. For samples irradiated by electrons or light, 
the Bourgoin mechanism described in previously can lead to 
athermal diffusion. On the other hand, radiation damage from 
ions can deposit enough energy in a small volume to directly 
create clusters of defects or even amorphous pockets and local 
melting despite cryogenic background temperatures. 

Besides the possibility of clustering even at low temper- 
atures, V — V interactions may provide another possible ex- 
planation of the low observed Frenkel pair formation volume. 
The largest supercell in our calculations that stabilized the T c j 
symmetry was the 128-site supercell, having //[V^] = 3.51 eV 
and v f [V^] = 0.34vs;, (Figure 4), with larger values of each 
expected with a larger supercell size. This would imply 
that h f [V +I} = 7.21 eV and v f [V + 1] = -0.24v Si . Inter- 
estingly, summing with jjjm) instead of Ij + would lead to 

h f [V +I} = 7.26 eV and v f [V +1} = -0.04v si . A single defect 
in a 100-site supercell corresponds to a defect concentration 
of 1%. At sufficiently high vacancy concentrations, V — V 
elastic interactions may destabilize the Jahn-Teller distortion, 
leading to a positive vacancy formation volume and a small 
Frenkel pair formation volume, as well as a lower formation 
enthalpy. The critical concentration to destabilize the Jahn- 
Teller distortion for actual vacancies in silicon need not lie in 
the interval jig < cy < jYg • 

Some earlier experiments claimed to find that cyDy = cjDj 
at T = 800 °C [5] or T = 1000°C [2]. Near that temper- 
ature, the effective activation enthalpy of Si self-diffusion 
should be close to 5(/i d [V] +/i d [/]). Our calculations indicate 
h f [V] = 3.69 eV and h f [I] = 3.68 eV under conditions of no 
doping or applied pressure. Our estimate of h m [V] = 0.27 eV 
appears reasonable, but we take h m [I] = 0.22 eV based on a 
study of low-symmetry pathways not investigated here [21]. 
Thus we arrive at h A [V] = 3.96 eV and h d [I] = 3.90 eV. These 
estimates are rather lower than most experimental reports. 
DFT methods predict migration enthalpies more accurately 
than they do formation enthalpies. Probably the most ac- 
curate method applied to the calculation of intrinsic defects 
in Si is diffusion quantum Monte Carlo (DMC), which gives 
h f [I%] = 4.82 eV [22], 1 .02 eV greater than the PW91 value, 
and leads to h d [I] =5.04 eV. 



PRESSURE-FERMI LEVEL STABILITY DIAGRAMS 

Examining the combined effects of changing pressure and 
Fermi level requires finding the enthalpies and volumes of for- 
mation of electrons and holes in Si as a function of pressure. 
The enthalpies are essentially the conduction band minimum 

hl(P)=H[Si^}(P)-Nh Si (P) (3) 
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FIG. 6: Minimum formation enthalpy of defects as a function of 

P and %). The intrinsic Fermi level is indicated as a gray line. 



and the valence band maximum 

h{(P)=Nh Si (P)-H[Si+ 1 ](P). (4) 
and the volumes are 

vl(P) = V[Si N l ](P) -Nv Si (P) = (5) 

and 

-v{(P) =N Vsi (P)-V[Si+ l ](P) = d -^p-. (6) 

At P = 0, we have v f e = 0.679v S i and v{ = -0.788 v S i, a siz- 
able effect. 

Since we have the pressure dependence of all these quanti- 
ties, we can essentially sweep Figure 5 over a finite pressure 
range, where the domain of the plot is the band gap. Showing 
all the enthalpy surfaces obscures the most important informa- 
tion, the lowest enthalpy defects. More useful is a plot of just 
the minimum enthalpy required to form a defect at a particu- 
lar pressure and Fermi level (Figure 6). The vacancy surface 
(Figure 6(a)) shows clear changes with pressure and marked 
curvature, indicating a variable formation volume, while the 
shape of the interstitial surface (Figure 6(b)) is close to being 
a prism. The location of the Fermi level is important for in- 
terstitials, but the pressure is not. Projecting these surfaces 
down onto the (P,£p) plane reveals other differences. The 
vacancy stability diagram (Figure 7(a)) displays a wealth of 
features: the instability of + under pressure, the transition 
from V£ to Vg , and in general the large changes in vacancy 
levels caused by pressure. None of these features are evident 
for interstitials (Figure 7(b)). 
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FIG. 7: P — Bp stability diagram of intrinsic defects in Si formed 
by projecting the surface of minimum enthalpy in Figure 6 for Vl 
(red), Vg (green), I T (red), I H (green), and /(no) (blue). The intrinsic 
Fermi level is indicated as a gray line. 

RELAXATION VOLUME TENSORS 

The stress state in integrated circuits is seldom hydrostatic. 
On the scale of the wafer it is typically biaxial due to surface 
oxidation or nitridation, but at a sub-micron scale becomes 
quite complicated, with enormous stress gradients. Aziz [23] 
has described the effect of the formation volume of defects in a 
film under biaxial stress, but points out that a full treatment re- 
quires more information about a defect than just its scalar vol- 
ume. Many calculations of defects in silicon point out changes 
in the positions of the neighboring atoms, but this can provide 
only qualitative, not quantitative information about the com- 
plete elastic distortion caused by the defect. The present work 
fills this gap. 

Among the features of the code used here is full relaxation 
of the vectors defining the supercell. This tells us how the 
shape as well as the size of the supercell change when a defect 
is formed. Were we dealing with a simple cubic cell changing 
isotropically from edge length Lq to L, we could define the 
relaxation volume as 

Instead, we have three non-orthogonal vectors (forming a 
transformation matrix L = Ly) changing in arbitrary direc- 
tions. We may generalize (7) as 

v rel =det(L )ln(L 1 L) ^de^Lo)^ 1 (L-Lq 1 ) (8) 

The relaxation volume tensors defined in this way are 
equivalent to the strain dipole tensors 

Ay = v#/v si (9) 
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TABLE IV: Relaxation volume tensors Vy of neutral defects in 

Si at Gjj = and £p = £;. The scalar part is indicated as well as 
the eigenvalues and corresponding principal directions. Scalar vol- 
umes in parentheses reflect the relaxation volumes calculated using 
the boundary conditions described previously (i.e. purely dilatational 
strain of the supercell). 



vacancy 


ViO'si) 


n, 


interstitial 


V{( v Si) 


n, 


v° L 

-1.02 
(-1.07) 


-0.83 
-0.82 
+0.63 


[110] 
[110] 
[001] 


1° 

(HO) 
+0.68 

(+0.59) 


+0.57 
+0.09 
+0.02 


[110] 
[001] 
[110] 


v° B 

-1.30 
(-1.30) 


-0.95 
-0.18 
-0.18 


[111] 

[211] 
[011] 


7° 
'h 

+0.64 
(+0.62) 


+0.23 
+0.23 
+0.19 


[2TT] 
[011] 

[111] 



4> +0.12 [111] 

+0.37 +0.12 [2TT] 

(+0.37) +0.12 [011] 



symmetry, drawing all six of its neighbors in, mostly toward 
each other. 

Among the self-interstitials, the Ij is isotropic, and the ifj 
is nearly so, pushing out a bit more in the (111) plane. The 
greatest component of the displacement around ^ 110 ) i s > as 
expected, an expansion in the [110] direction, but the other 
components are quite small. Two such interstitials, aligned in 
the same direction, will tend to repel if their axes are parallel 
to the line between them, and attract if the axes are perpendic- 
ular to the line between them. This suggests that clusters of 
these interstitials may agglomerate in a {110} plane with their 
axes aligned to maximize their attraction. This long-range in- 
teraction may explain how interstitials are drawn together to 
form extended defects, including (110) chains and eventually 
{311} defects, as seems to be the case [27]. 



PREVIOUS WORK 



discussed by Nowick and Berry [24] and closely related to the 
piezo-spectroscopic elastic dipole tensors 

Pij = C ijkl v£l (10) 

of Kroner [25] (where Cy# is the elastic modulus tensor of the 
material). 

Using (8), we have calculated the full tensor volumes of 
relaxation of intrinsic defects in Si, which are mostly quite 
anisotropic. The results are more easily interpreted by diag- 
onalizing the matrices to find the eigenvalues and principal 
directions, which we have compiled in Table IV. (Some of 
the defects have an axis of symmetry, leading to degeneracy 
of eigenvalues. As a result, the corresponding eigenvectors 
are not unique, and a different pair of directions in that plane 
could be chosen instead.) Useful methods of visualizing these 
tensors [26] include plotting the volume ellipsoids (Figure 8) 
defined by 

x 2 y 2 z 2 

- + - + - = 1 (ID 

vf v 2 v 2 

and the volume director surfaces (Figure 9) 

2 2 2 

X -+ y - + Z - = ±l (12) 

VI V2 V 3 

where +1 is taken for expansion and —1 for contraction. 

The sum of the eigenvalues — the scalar part of the tensor- 
is approximately equal to the scalar relaxation volume calcu- 
lated previously. The principal directions and the signs of the 
eigenvalues are all in accord with intuition. We find that the 
neutral lattice vacancy has Z??^ symmetry (at least at this pres- 
sure). The expansion in the [001] direction makes its slightly 
negative scalar relaxation volume all the more noteworthy. 
The degree of anisotropy makes clear that scalar relaxation 
volumes are far from a complete picture of the defect dis- 
placement field: The split vacancy has a quite different Ct, v 



First-principles calculation of the effect on intrinsic defect 
formation energies of changing lattice parameter was per- 
formed as early as 1984 by Car et al. [28] and LDA in 1989 by 
Antonelli and Bernholc [29], using supercells with 32 atoms. 
Those works did not include /(no)- Sugino and Oshiyama 
[30] studied the effect of pressure on the diffusion of group 
V dopants (P, As, Sb) in silicon, but not silicon self-diffusion. 
Most published reports are limited to examining energy dif- 
ferences with a fixed volume. 

Previous work by Zhu [31] and others [32] used the local 
density approximation (LDA), which does not include the gra- 
dient correction of GGA (generalized gradient approximation) 
methods such as PW9 1 . GGA methods tend to predict ener- 
gies of localized states (such as defects) with less error than 
LDA. Also, those calculations were limited to 64-site cells. 
However, it has been found necessary to use supercells con- 
taining over 200 atoms to stabilize the Jahn-Teller distortion 
ofVf (cf. Puskaefa/. [15]). 



CONCLUSIONS 

We have performed first-principles calculations on a num- 
ber of basic properties of intrinsic defects in silicon, some 
novel, some not presented together in a unified analysis be- 
fore. The relaxation volumes of electrons and holes (+0.68vsi 
and — 0.79vsi, respectively) are an appreciable fraction of 
an atomic volume in magnitude. The formation enthalpy 
of a neutral vacancy is 3.69 eV and its migration enthalpy 
is 0.27 eV. The relaxation volume of the neutral vacancy, 
— 1.07vsi, is of the expected sign but the magnitude of 
the number is rather large, resulting in a formation vol- 
ume of — 0.07vsi, with an activation volume of migration of 
— 0.24v's; . That is, hydrostatic pressure should lead to a slight 
increase in equilibrium vacancy concentration and an increase 
in vacancy diffusion. The most stable self-interstitial species, 
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FIG. 8: Volume ellipsoids of neutral defects in Si at cry = and £f = £;, as defined in Eq (1 1). Principal axes are indicated in gray. 




FIG. 9: Volume director surfaces of neutral defects in Si at <jy = and £p = £j, as defined in Eq (12). Surfaces of both expansion (light 
blue) and contraction (red) are shown. Principal axes are indicated in gray. 



Ij , has a formation enthalpy 3.68 eV, equal to the forma- 
tion enthalpy of V® within the accuracy of our method, and a 
formation volume of — 0.57vs; (v rel = +0.43vs;). 

To estimate the concentration of a defect requires calcula- 
tion of its vibrational and configurational entropy. The calcu- 
lations we have performed do not provide this information, but 
the near equality of the formation enthalpies of the most sta- 
ble vacancy and self-interstitial demand that the ratio C//Cy 
of their equilibrium concentrations should not vary much with 
temperature. The ratio will depend upon pressure, however. 
The formation volume of a Frenkel pair with no clustering or 
amorphization should be a sizable negative number, — 0.48vs; 
to — 0.68vsi, unless interactions at high vacancy concentra- 
tions destabilize the Jahn-Teller distortion. These parameters 
should be of interest in studies of silicon under large elastic 
stresses. They imply that increasing hydrostatic pressure in- 
creases the equilibrium concentration of both vacancies and 
self-interstitials, though the effect on vacancy concentration 
should be weak. Increasing pressure should also increase the 
mobility of vacancies. 

We have presented stability diagrams (akin to phase dia- 
grams) of the intrinsic defects in silicon, showing the com- 
plexity of vacancies under pressure, in stark contrast to the 
behavior of self-interstitials. We have verified Bourgoin's pre- 
diction of a crossing in the vacancy potential energy surfaces. 



Finally, we have calculated the full tensor relaxation vol- 
umes of these intrinsic defects, enabling researchers to model 
the biased diffusion of both vacancies and self-interstitials un- 
der the non-hydrostatic stress states found in actual devices. 

This work was performed under the auspices of the 
U. S. Department of Energy by the University of California, 
Lawrence Livermore National Laboratory under Contract No. 
W-7405-Eng-48. Funding was provided by the Department of 
Energy Office of Basic Energy Sciences. 
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